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Introduction 


This  work  focuses  on  constructing  and  solving  the  initial 
value  problem  for  systems  of  fractional  order  differential 
equations.  These  equations  arise  when  fractional  order 
derivatives  are  used  to  describe  the  linear  viscoelastic  behavior 
of  elastomeric  materials  in  the  equations  of  motion  of  damped 
structures  (1) (3) (5) .  The  fractional  order  derivative  is 
particularly  well  suited  to  describe  the  frequency  dependent 
stiffness  observed  in  many  viscoelastic  materials  (2) (6) 
( 7 ) ( 8) ( 9 ) ( 20)  .  These  viscoelastic  stress  -  strain  constitutive 
relationships  are  constructed  by  replacing  the  ordinary 
derivatives  appearing  in  the  classical  linear  viscoelastic  model 
(11)  with  fractional  order  derivatives  (6) (8) (14) (15) (18)  . 

Because  of  their  ability  to  describe  weak  frequency 
dependence,  significantly  fewer  fractional  derivative  terms  are 
needed  to  relate  time -dependent  stress  and  strain  fields. 
Typically  models  with  four  or  five  parameters  effectively  describe 
the  material  stiffness  over  five  or  more  decades  of  frequency. 
The  small  number  of  parameters  makes  least  squares  fitting  of  the 
data  attractive  and  several  materials  have  been  modeled 
(3) (4) (6) (20) .  These  constitutive  models  are  empirical,  but  their 
general  mathematical  form  is  suggested  by  approximate  molecular 
theories  that  relate  microstructure  to  macro  viscoelastic 
stiffness  properties  (4) . 

More  important  to  the  engineer,  these  models  can  be 
straightforwardly  incorporated  into  continuum  and  finite  element 
formulations  for  the  analysis  of  structures  containing  in 
principle  and  unlimited  number  of  linear  elastic  and  linear 
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viscoelastic  materials  (1).  These  formulations  lead  to  closed 
form,  real,  continuous  and  causal  solutions  using  Laplace 
transforms  or  Green's  function  solutions  and  convolution  (1)(3). 
The  fractional  calculus  model  that  produces  these  solutions,  like 
most  viscoelastic  models,  is  a  hereditary  operator  that  requires  a 
complete  displacement  or  strain  history  to  prepuce  precise 
results.  Present  solution  techniques  are  predicated  on  a  total 
history  (1) (3) (5) . 

Fortunately,  the  fractional  calculus  model  is  also  a  fading 
memory  operator  where  events  in  the  distant  past  have  less  effect 
on  the  present  and  future  states  than  do  comparable  events  in  the 
recent  past.  This  feature  raises  the  possibility  of  ignoring 
events  in  the  distant  past  and  generating  approximate  structural 
responses  based  the  recent  past  (a  few  characteristic  times  back 
in  history),  the  present  state  and  future  loadings.  When  the 
effects  of  the  recent  past  are  cast  as  pseudo  forcing  functions 
superimposed  on  the  future  loads,  one  has  constructed  an  initial 
value  problem  where  the  structural  response  is  uniquely  dependent 
on  the  present  state:  and  future  loads.  Solving  the  initial  value 
problem,  instead  of  calculating  the  response  based  on  a  complete 
displacement  history,  clearly  reduces  the  effort  required  to 
produce  solutions  and  obviates  the  need  to  begin  the  analysis  in 
the  distant  past  at  which  time  the  viscoelastic  materials  were  in 
a  virgin,  quiescent  state. 

The  initial  value  problem  is  posed  as  a  system  of 
integro-differential  equations,  which  upon  close  examination  may 
be  viewed  as  differential  “quations  generalized  to  fractional 
order.  These  equations  are  referred  to  as  fractional  order 
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differencial  equations  with  some  justif ication  because  of  their 
similarity  with  ordinary  differential  equations  with  constant 
coefficients.  For  instance,  it  is  shown  that  a  system  of  m 
fractional  order  differential  equations  produces  m  eigenfunctions 
needing  m  initial  conditions  for  a  unique  solution.  These 
eigenfunctions  will  be  cast  in  terms  of  Mittag-Lef fler  functions 
(16),  long  viewed  by  some  as  generalized  exponential  functions 
(12 : 280) (15:260) (18:44) .  The  total  solution  for  a  fractional 
order  differential  equation  is  seen  to  be  composed  of  a 
generalized  homogeneous  solution  uniquely  dependent  on  the  initial 
value  and  a  generalized  particular  solution  uniquely  dependent  on 
the  forcing  function. 

Setting  the  fractional  order  to  integer  values  produces 
ordinary  differential  equations  with  constant  coefficients 
(13:527)  and  the  generalized  solutions  become  the  traditional 
solutions  to  systems  of  ordinary  differential  equations.  In 
addition  it  is  shown  that  the  fractional  order  differential 
equations  produce  non-singular  Green's  functions  which  insure 
continuous  dependence  on  the  data  which,  when  coupled  with 
existence  and  uniqueness  considerations,  lead  to  well-posed 
problems.  The  preponderence  of  the  strong  parallels  between 
fractional  order  differential  equations  and  ordinary  differential 
equations  lends  credence  to  the  view  that  this  practical  tool  for 
the  engineer  is  in  fact  a  generalization  to  fractional  order  of 
the  theory  of  ordinary,  linear,  differential  equations  with 
constant  coefficients. 
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Constructing  the  Fractional  Order  Differential  Equations  of  Motion 
Differential  equations  of  fractional  order  arise  in  solid 
mechanics  when  rational  fractional  order  derivatives  are  employed 
to  describe  the  viscoelastic  phenomenon.  This  approach  produces 
fractional  order  differential  equations  of  motion  where  the 
elastic  constants  are  replaced  by  fractional  derivative  operators 
intended  to  describe  the  time  -  dependent  viscoelastic  stress - s train 
moduli.  The  extended  Riemann  Liouville  fractional  der ivative ^ 

( 17 : 50) ( 19 : 19 )  is  a  linear  operator  and  is  defined  as 


V)[u(t) 


,t  - H(t) -  dr 

Jo  r(l-a) (t-r)Q 


0  <  a  <  1 


This  operator  is  used  to  construct  viscoelastic  stress  -  strain 


constitutive  relationships  of  the  form 


c'(t,xi) 


E  <r(t,x.)  +  S  E  D/.  [« (t,x.)] 

o  ’  i'  p  (t) ‘  i/ 1 


shown  here  for  uniaxial  deformation.  This  constitutive 
relationship  may  be  viewed  as  a  generalization  to  fractional  order 
(1:17)  of  the  classical  viscoelastic  model  (11:14)  where 
derivatives  of  integer  order  are  used  to  relate  empirically 
time-dependent  stress  an!  strain  fields. 

This  fractional  derivative  constitutive  relationship  has  been 
successfully  incorporated  into  the  equations  of  motion  for 
continuum  and  finite  element  formulations  producing  real, 


This  definition  is  usually  presented  with  a  restricted  to  less 
than  one.  However,  Giittinger  (13:527)  has  proven  that  the  weak 
limit  of  the  kernel  as  a  tends  to  one  is  the  Dirac  delta  function, 
in  which  case  the  definition  produces  a  first  derivative. 
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continuous  and  causal  solutions  for  forced  structural  response 
The  solution  method  for  the  resulting  equations  of  motion  is  based 
on  decomposing  them  into  systems  of  decoupled  linear,  fractional 
order  differential  equations  (3)(5). 

The  overall  approach  may  be  demonstrated  by  constructing  the 
equation  for  forced  motion  of  a  simple  spring-mass  harmonic 
oscillator  and  generalizing  the  result  for  multiple  degree  of 
freedom  structures  having  viscoelastic  components. 

F(t)-Mu(t)=ku(t)  (3) 

Modeling  the  spring  as  a  massless  bar  (with  elastic  modulus  E, 
cross  -  sectional  area  A  and  length  L)  yields 

a(t)  -  F(t)  -  Mu(t)  -  Eu( t)  =  Ee(t)  (4) 

A  L 

This  may  be  viewed  as  an  elastic  stress  -  strain  relationship  where 
the  inertially  induced  stress  is  described  using  D'Alembert's 
principle.  The  massless  rod  is  now  taken  to  be  viscoelastic  and 
modeled  using  eqn  2,  where  the  elastic  terms  and  the  first 
fractional  derivative  terms  acting  on  stress  and  strain  are 
retained . 

(1  +  bDQ)a(t)  =  (Eq  +  E1D®)€(t)  (5) 

The  resulting  equation  of  motion  in  operator  form  is 


(1  +  bDQ) 


F( t )  -  M  P2u( 
A 


Using  the  composition 
derivative  operator, 


t) 


=  (E  +E,D  ) 
o  1 


u(t) 


property  (17:30) 


of  the 


(6) 

fractional 
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the  equation  of  motion  takes  the  form 


bMD2+au(t)  +  MD2u(t)  +  k1D°u(t)  +  k  u(t)  =  bDQF(t)  +  F(t).  (8) 

Here  the  equations  of  motion  have  been  expanded  to  describe 

the  behavior  of  systems  with  several  degrees  of  freedom.  The 

double  and  single  underlined  variables  are  square  and  column 

matrices,  respectively,  of  order  equal  to  the  number  of  physical 

degrees  of  freedom,  N.  Identifying  the  largest  factor  /3 ,  of 

the  form  1/n,  common  to  all  the  fractional  orders  of 

2 

differentiation  in  eqn  8  and  again  applying  the  composition 
property,  this  equation  of  motion  becomes 


0.m 


bM  ( D  ) 


m(dV 


ki(D^)q 


k 

“OJ 


\  o' 

u(t)  -  1  +  b(u  )q 

J 


F(t) 


(9; 


3  m 

where  m,  w  and  q  are  integers  and  (D;  is  the  0  order  derivative 
taken  m  times.  The  orders  of  the  corresponding  differential 
operators  in  eqns  8  and  9  are  also  equal . 

/3m  =  2  +  q 

-  2  (10) 

Pq  =  a 

0  =  1/n 

The  most  general  form  of  these  equations  of  motion  is 


This  restriction  on  0  is  tied  to  solving  the  initial  value 
problem.  When  solving  for  the  motion  starting  with  zero 
displacement  and  velocity,  this  restriction  is  not  needed  and  /3 
may  be  taken  to  be  the  Targest  rational  factor  common  to  all 
orders  of  differentiation. 
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(11) 


l  c  (dV  u ( t )  =  (1  +  b(D^)q)F(t)  =  f*(t). 

p=0  =P  ~ 

Here  the  c  are  real,  although  many  may  be  zero,  and  f  (t)  is  the 
=P 

result  of  the  viscoelastic  stress  operator  acting  on  the  applied 
forces,  F(t),  as  shown  in  eqn  9. 

Eqn  11  poses  the  system  as  N  equations  of  order  /9m  that  can 
be  alternatively  posed  as  m-N  equations  of  order  /?.  In  matrix 
form  the  m-N  equations  of  0  order  are 


'  0 

0 

A 

0 

C 

L  =m 

: — i 

H 

Oil 

CO 

Oil 

CO 

OH 

(D^)m_ 1u(t ) 


(DVu(t) 

(DVu(t) 


u(t) 


(12) 


0  ‘ 

(D^)m— 1u(t)' 

0  1 

-A 

• 

+ 

0 

«» 

(D^)Zu(t) 

►  =  - 

0 

0 

(DVu(t) 

0 

_  0 . . 0  0 

C 

=0J 

u(t) 

l  £  (t)J 

where  the  lowest  equation  is  seen  to  be  eqn  11.  The  matrix  [A]  is 
chosen  such  that  both  square  matrices  of  order  m-N  become 
symmetric  and  the  top  (m-l)-N  equations  are  satisfied  identically. 
This  is  accomplished  by  constructing  A  such  that  all  matrices,  c^ 
lying  on  any  given  diagonal  running  from  lower  left  to  upper  right 
in  the  first  matrix  of  eqn  12,  are  equal.  This  form  of  the 
equations  of  motion  will  be  referred  to  as  the  expanded  equations 
of  motion. 

For  example,  if  a  is  one  half  in  eqn  8,  then  0  is  one  half 
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making  m-5  in  eqn  9 ,  and  the  expanded  equations  of  motion  become 


'  0  0  0  0  bM  ' 

'(DWVu(t)' 

0  0  0  bM  M 

(DW2)3u(t) 

0  0  bM  M  0 

* 

(D1/2)2u(t)- 

0  bM  M  0  0 

(D1/2)1u(t) 

bM  M  0  0  k 

u(t) 

_  _  _  _  -i 

> 

0  0  0  -bM  0 

'(D1' Vu(t)' 

f  \ 

o 

0  0  -bM  -M  0 

(D1/2)3u(t) 

o 

+ 

0  -bM  -M  00 

» 

(DW2)2u(t) 

-  =  •. 

0 

-bM  -M  0  0  0 

(DW2)1u(t) 

0 

o 

^11 

oil 

oil 

oil 

oil 

u(t) 

i(l+bD1/2)F(t) 

(13) 


* 

f  (t) 


Both  the  general  form,  (eqn  12),  and  the  example  in  eqn  13 


are  now  posed  in  terms  of  two  real,  square,  symmetric  matrices  for 


which  an  orthonormal  transformation  exists 


\ 

(D^)2u(t) 

$ 

- 

y3(t> 

(D/?)1u(t) 

L 

y2(c) 

u(t) 

y 

ypo 

/ 

that  leads  to  a  system  of  m-N  uncoupled  differential 


order  p. 


D 


$ 


equations  of 


(15) 


Moreover,  one  can  see  that  this  process  is  analogous  to 
posing  systems  of  higher  older  coupled  ordinary  differential 
equations  with  constant  coefficients  as  a  collection  of  first 
oi'der  differential  equations,  the  difference  being  that  fractional 
order  differential  equations  result  here. 
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At  this  Doint  one  might  ask  why  we  have  resorted  to  this 
formalism.  In  general,  the  stiffness  matrices  and  kQ  are 

distinctly  different  matrices  when  two  or  more  materials  are 
present  in  the  structure  and  at  least  one  is  viscoelastic. 
Constructing  an  orthonormal  transformation  that  decouples  a  system 
with  more  than  two  (in  this  case  two  stiffness  and  one  mass) 
symmetric  matrices  is  generally  not  possible.  Consequently, 
casting  this  system  of  equations  in  terms  of  two  symmetric 
matrices,  eqn  12,  becomes  necessary  to  construct  manageable 
decoupled  equations  through  an  orthonormal  transformation. 

The  major  drawback  of  having  to  rely  on  the  formalisms  of 
eqns  12  and  14  is  that  they  require  the  manipulation  of  large, 
usually  sparsely  populated,  matrices.  Numerical  methods,  that 
capitalize  on  the  repetitive  structure  of  eqn  12  and  the  less 
obvious  repetitive  structure  of  eqn  14  to  increase  efficiency, 
have  been  applied  to  extract  eigenvalues  and  eigenvectors  for  a  65 
degree  of  freedom  system  with  only  modest  computational  effort 
(5:923).  When  using  fractional  derivative  models  to  determine  the 
response  of  structures  containing  viscoelastic  components,  the 
formalism  of  eqn  12  is  necessary,  if  decoupled  equations  of  motion 
are  required  for  modal  analyses.  In  any  event,  the  fractional 
order  basis  equations  shown  in  eqn  15  are  the  progenitor  of  the 
initial  value  problem. 
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Constructing  the  Initial  Value  Problem 

The  decoupled  fractional  order  equations  of  motion  or  basis 
equations,  (eqn  15)  individually  Lake  the  form 

(D^+  a^)y(t)  =  f(t)  P  =  1/n  (16) 

where  the  subscripts  have  been  dropped  to  simplify  notation. 
Green's  function  solutions  for  these  equations  are  relatively 
straightforward  and  the  resulting  expressions  for  the  forced 
response  of  the  structure  can  be  shown  to  be  real,  continuous  and 
causal  (1:73)  These  solutions  to  eqn  16  may  be  viewed  as 
particular  solutions  of  the  fractional  order  differential  equation 
of  motion. 

It  is  curious  to  note  that  the  only  homogeneous  solution  to 
eqn  16  is  the  trivial  solution.  This  appear?  to  be  consistent 
with  a  strict  interpretation  of  eqn  2,  the  generalized 
viscoelastic  constitutive  model.  Inherent  to  the  model  is  the 
implication  that  the  relationship  between  stress  and  strain  at 
any  given  time  should  be  a  function  of  the  complete  stress  and 
strain  histories.  In  other  words,  at  time  zero  the  viscoelastic 
material  should  be  in  its  virgin,  undeformed  state  where 
structural  motion  is  commencing  from  a  quiescent  state. 
Attempting  to  impose  non-trivial  initial  conditions  implies  the 
existence  of  previous  motion  that  is  inconsistent  with  the 
viscoelastic  model  and,  therefore,  homogeneous  solutions  are  not 
needed . 

In  practice  one  usually  does  not  have,  cannot  calculate  or 
should  not  calculate  a  complete  time  history  of  structural  motion 
to  conform  with  this  restriction.  For  example  consider  the  case 
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of  a  viscoelastically  damped  structure  subjected  to  several 
episodic  loadings,  where  time  intervals  between  episodes  do  not 
allow  the  structure  to  assume  an  essentially  undeformed,  quiescent 
state.  As  a  practical  matter,  one  would  like  to  be  able  to 
determine  the  structure's  response  for  one  loading  episode  where 
its  initial  displacements  and  velocities  are  significant  and 
known,  and  a  "recent"  history  of  structural  motion  is  available. 

lo  produce  equations  of  motion  which  satisfy  this  need,  eqn 
16  will  be  posed  in  terms  of  a  shifted  time  scale,  t,  shown  in 
Figure  1. 


1 


T(l-/3) 


dt 


y(r-t  ) 


o  (t-r) 


0  dr  +  a  y(t-to)  = 


(17) 


Here  t  is  the  time  of  the  onset  of  the  loading  history  for  which 
the  response  is  needed.  Time  zero  on  the  t  timescale  is  taken  to 
be  several  characteristic  times,  a  \  preceding  the  beginning  of 
this  load.  This  preceding  time  interval  is  not  restricted  to 
being  load  free  and  in  general  may  have  loading  present  right  up 
to  the  start  of  the  loading  history  of  interest.  The  loads  prior 
to  tQ  are  f  (t-tQ)  and  the  corresponding  response  is  y(t-t  ). 
The  equation  of  motion  resulting  from  f(t-t  )  is 

D  y( t— tQ)  +  a/3y(t-to)  =  f(t-tQ)  (18) 

The  loads  for  the  episode  of  interest  (t  >  t  )  are  ?(t-t  )  and 

o  ~ 

the  equation  for  the  corresponding  response  y(t-t  )  is 

D^y ( t- t  )  +  a^y(t-t  )  =  ?(t-t  ) 
o  o  o 

The  total  response  for  t  >  t  is  y  +  y  and  the  general 


1* 


(19) 


expression  for  the  response  is 


r(l-/9) 


y* (r-u)  +  y'  (r-u)  du  +  aP 


y(r)  +  y(r) 

=  ?(r)  +  g(r) 


(20) 


where  r  —  t-t  ,  U  =  t’r‘  Here  g(r)  is  a  pseudo  forcing  function 
that  produces  the  residual  response  of  the  structure  due  to  the 
prior  application  of  f(t-t  ). 


g(r) 


r(i-/9) 


r+  t 

o 

y'  (r-u) 


u 


p 


du  + 


y(-V 


(r+t  ) 

o 


/9 


(21) 


Expressing  eqn  20  in  terms  of  the  t  time  scale  in  Figure  1,  where 
zero  time  is  now  the  onset  of  the  loading  episode  of  interest, 
yields 


r(l-/0) 


y  J.5 ~T)dr  +  a^y(t)  =  ?(t)  +  g(t)  -  g(t). 

0  TP 


(22) 


Note  that  here  the  order  of  differentiation  and  integration 
in  the  fractional  derivative  operator  is  the  opposite  of  eqn  1. 
This  reversal  of  operations  occurred  when  Leibnitz's  rule  was  used 
to  differentiate  the  integral  in  eqn  17,  producing  eqns  20  and  21. 
This  change  will  prove  crucial  to  solving  the  initial  value 
problem,  because  in  contrast  with  eqn  16,  eqn  22  possesses  both  a 
particular  solution,  uniquely  dependent  on  the  forcing  function, 
and  a  homogeneous  solution,  uniquely  dependent  on  the  initial 
value ,  y (0) . 

Before  presenting  these  solutions  it  is  important  to  address 
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Figure  1  -  Time  Scales  for  the  Loads  and 
Responses  of  the  Initial  Value  Problem. 


the  relationship  between  the  operator  appearing  in  eqn  22  and  the 
original  definition  shown  in  eqn  1.  Using  Leibnitz's  rule  to 
differentiate  the  integral  in  eqn  1  yields 


1 

_d  rL 

u(t'T)  dr 

1  1 

I  u(0) 

-t 

u 

(t'°  dr1 

L  (  ?  7  'i 

r(i- 

a) 

dt  „ 

a 

r( 

1-a)  j 

1 

a  dr 

r  J  / 

0 

T 

J  0 

7  J 

or  in 

oper 

ator 

form 

Da‘ 

u(t') 

u(0)t'Q 

r(i-o) 

A 

+  dq 

u(t)] 

(24) 

where 

DG 

is  the  defini 

tion 

A 

and  D  is 

the 

modified 

operator 

appear 

ing 

in 

eqn  22. 

In 

fact 

D°  is 

the 

Riemann- 

Liouville 

fractional  integral  of  order  1-a  of  the  first  derivative  of  the 

3 

function  or  effectively  an  order  -a  integal  of  a  function.  In 
effect  the  operator  DQ  treats  u(t)  as  though  it  is  zero  for 
negative  time  and  is  "turned  on"  by  a  step  function  at  time  zero. 
The  singular  term  appearing  in  eqns  23  and  24  is  the  a  order  time 
derivative  of  a  step  function  with  magnitude  u(0) .  Conversely, 

A 

Cl  + 

the  D  operator  treats  u(t)  as  though  its  value  for  t  =  0  is  an 
analytic  continuation  of  its  non-zero  value  at  t  =  0  .  The 
existence  of  an  analytic  continuation  of  u(t)  into  positive  from 
negative  time  means  its  first  derivative  is  bounded  at  time  zero. 
As  shown  in  Appendix  C,  a  bounded  first  derivative  at  t  -  0 
restricts  the  initial  values  of  the  modified  fractional 
derivatives  of  positive,  rational  order  less  than  one  to  zero. 
This  result  will  prove  crucial  to  the  solution  of  the  initial 


The  relationship  between  the  operators  (eqn  24)  is  developed 
more  formally  by  Oldham .( 17 : 50) 
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value  problem. 

Posing  equation  22  in  terms  of  the  modified  fractional 
derivative  operator,  D  , 

(D^  +  a?)  y ( t)  =  ?(t)  +  i(t)  =  g(t)  (25) 

produces  the  modified  basis  equations.  Note  the  similar 

appearance  of  eqns  16  and  25.  Recall  that  eqn  16  is  based  on  the 
t  time  scale  and  has  a  trivial  homogeneous  solution.  On  the  other 
hand,  eqn  25  is  based  on  the  t  time  scale,  possesses  a  non-trivial 
homogeneous  solution  and  accounts  for  the  effects  of  previous 
motion  through  the  initial  value,  y(0) ,  and  pseudo  forcing 
function,  g( t) . 
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Solving  the  Initial  Value  Problem 


The  eventual  goal  is  to  use  the  solutions  of  the  modified 
basis  equations,  eqn  25,  to  construct  solutions  to  the  original 
structural  equations  of  motion  for  non-trivial  initial 
displacements  and  velocities,  where  the  relaxation  effects  induced 
by  previous  motion  are  accounted  for  by  the  pseudo  -  fore ing 
functions,  as  in  eqn  20.  The  overall  homogeneous  solution  will  be 
a  superposition  of  the  homogeneous  solutions  for  the  modified 
basis  equations  and  will  be  shown  to  satisfy  the  initial 
conditions.  The  overall  particular  solution  will  be  constructed 
from  the  particular  solutions  to  the  modified  basis  equations, 
derived  using  Green's  functions.  Superimposing  the  overall 
homogeneous  and  particular  solutions  produces  the  total  solution 
to  the  initial  value  problem  for  structural  motion. 

The  overall  homogeneous  solution  is  constructed  by  first 
solving  eqn  25  for  all  the  homogeneous  solutions  to  the  modified 
basis  equations.  These  solutions  take  the  form 


Vc) 


n,<°>  1 


t 


(at)  V 


P-0 


r(i+p/3) 


(26) 


which  is  a  special  case  of  the  beta  order  Mittag-Leffler  function 
defined  as  (16 . 102) 


£  (x)  =  y  _ (Jil 

e  PL  r(1 ' 

In  Mittag-Leffler  notation  the  homogeneous  solution  is 


(27) 


yh(t)  -  yh(0)  E0 


-(at)"  , 


(28) 
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where  this  special  form  of  the  Mittag-Leffler  function  has  the 
property 


D^E 


n 


-(at/ 


-a^E 


-(at)' 


(29) 


Including  the  particular  solution,  the  total  solution  to  each  of 
the  modified  basis  equation  is 


>’  ( t )  =  yh(0) 


g  (  t - t ) dr 


(30) 


which  can  be  determined  using  Laplace  transforms  or  other 
traditional  solution  techniques  for  5 ntegral - differential 
equations.  The  kernel  in  the  convolution  integral  of  eqn  29  is 
the  unit  impulse  solution  (Green's  function)  for  the  modified 
basis  equations,  and  is  singular.  Note  that  E^(0)  is  not  zero  and 
that  the  singular  behavior  of  the  kernel  can  be  determined  through 
a  straightforward  application  of  eqn  1. 

It  is  the  singular  nature  of  fractional  order  derivatives  of 
E^(-(at)  )  that  is  useful  in  resolving  an  apparent  paradox  in  the 
overall  initial  value  problem.  Recall  that  there  are  m-N  (eqn 
25)  modified  basis  equations  needed  to  characterize  the  structure, 
where  the  solution  for  each  modified  basis  equation  has  a 
homogeneous  solution  containing  a  different  initial  value.  For 
example  in  the  single  degree  of  freedom  problem,  N=1 ,  we  have  only 
two  initial  conditions,  displacement  and  velocity.  It  would 
appear  as  though  we  have  insufficient  information  to  combine  m 
homogeneous  solutions  to  the  m  modified  basis  equations  in  a 
unique  overall  homogeneous  solution. 

This  paradox  becomes  more  apparent  when  eqn  19  is  used  to 
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solve  for  the  m-N  initial  values  of  the  homogeneous  basis 
functions  in  terms  of  the  structure's  initial  displacements  u^(t) 
and  their  derivatives  evaluated  at  time  zero. 


f  \  ( 


(dV  1uh(t) 

. 

>m;N( 

(D^)2uh(t) 

•  = 

* 

y3(c) 

(D^) 1Hh(t) 

. 

y2(t) 

yx(t) 

The  paradox  is  that  at  this  point  only  u^(0)  arid  D  u^(0)  can  be 
specified,  while  the  remaining  elements  in  the  state  vector  on  the 
left  of  eqn  31  are  undetermined. 

It  is  curious  that  the  initial  value  problem  should  also  call 
for  the  initial  values  of  accelerations  and  higher  order 
derivatives  as  well.  Note  that  the  order  of  the  differential 
equations  of  motion  (eqn  8)  is  order  2+o  or  equivalently  /3m  and 
that  the  state  vector  in  eqn  31  calls  for  the  initial  values  of 
derivatives  up  through  2+a-/3  or  equivalently  /3(m-l).  In  other 
words,  when  posing  N,  /3m  order  differential  equations  as  a  system 
of  m-N  differential  equations  of  order  /3  the  corresponding  initial 
value  problem  calls  for  all  the  initial  values  of  the  p/3  order 
derivatives  of  the  displacement  vector,  u^(t) :  p  =  0,1,2,---,  m-1. 
These  requirements  appear  to  be  analogous  to  the  traditional 
initial  value  problem  for  ordinary  differential  equations,  but 
also  leaves  one  with  the  requirement  for  yet  more  initial 
conditions . 

It  is  proven  in  Appendix  C  that  all  of  the  non- integer 
derivatives  of  cl  (t)  of  order  less  than  two  appearing  in  the  state 
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vector  have  zero  initial  value.  The  initial  values  for 

acceleration  and  the  accompanying  higher  order  derivations 
appearing  in  the  state  vector  can  be  determined  by  returning  to 
the  original  equation  of  motion,  eqn  8,  and  using  successive 
applications,  of  eqn  24  to  determine  the  singular  terms  in  the 
equation  of  motion.  The  resulting  equation  of  motion  for  the 
response  to  turning  off  the  previous  forcing  function  is 


-  bM 


u(0)t 

To^) 


m- 2n- 1 
-  bM  l 
i-1 


(m-2n-f  )/3 

A 

D  u(O') 


+  (1  +  bD  )  M  u ( t ) 


u(°')t  “ 

Si  Fo^r  +  (^o +  ^iD 


-  b 


F(0  )t 

_r(l-a) 


+  G(t) 


(32) 


The  fractional  derivatives  in  this  equation  of  motion  are 
evaluated  for  t<0  or  equivalently  t<t  as  shovm  in  Appendix  A. 
G(t)  are  the  pseudo  forces  needed  to  produce  the  residual  motion 
associated  with  the  previous  loading  history,  already  accounted 
for  in  the  modified  basis  equations.  The  singular  forcing 
function  is  the  result  of  the  a  order  derivative  of  the  step 
funtion  turning  off  F(t) .  The  remaining  singular  behavior  is  the 
result  of  repeatedly  applying  eqn  24  to  separate  out  the  singular 
behavior  of  the  fractional  derivatives  of  acceleration. 

The  corresponding  equation  of  motion  for  the  response  to 
the  new  loads  is 
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u(0+)t'° 


bM 


=  r(l-a) 


m-  2n- 1  -i/3  ^(n-2n-£)fi 

+  bM  I  ■  —0,I>  3(0  ) 


i-1 


r(i-i/3) 


u(0+)t'Q 


(33) 


+  (1  +  bD  )M  u  ( t )  +  k  — =-tt — r  +  (k  +  k  D  )u(t) 
N  =  —  =1  r ( 1 - a )  =o  =1 


bF(0+)t"° 

Td-tt) 


+  (1  +  bDQ)F(t) 


where  the  singular  forcing  function  results  from  again  using  eqn 
2U  to  express  the  effects  of  the  step  function  turning  cn  F(t). 
The  remaining  singular  behavior  ir  also  the  result  of  using  eqn  24 
to  separate  out  the  singular  behavior  of  the  fractional 
derivatives  of  acceleration.  Again  the  tilde  and  double  tilde 
notation  identify  motion  due  to  previous  forces,  F(t),  and  present 
forces,  P(t),  respectfully,  as  in  eqns  18  and  19. 

Equating  the  coefficients  of  the  strongest  singularities 
(order  a)  in  eqns  32  and  33  yields  two  equations  needed  to 
establish  the  initial  conditions  on  acceleration. 


-  bM  u(O')  -  k1u(0')  =  -  bF(O')  (34) 

bM  u(0+)  +  k1  u(0+)  ~  bf(0+)  (35) 


Adding  these  two  equations  produces  the  relationship  needed  to 
establish  changes  in  the  initial  conditions  due  to  stopping  and 
starting  of  the  load  histories. 


u(0  )  -  u(0  ) 


+  b " 1k1  3(0+)  -  u(0  )|  -  F(0  )  -  F ( 0  )  (36) 


Since  this  relationship  is  based  on  step  loading,  which  is 
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incapable  of  instantaneously  changing  the  displacement  or  velocity 
time  history  between  time  0  and  0+ ,  one  can  conduce  that 

u(0+)  =  u(0_) 

u(0+)  -  u(O') 

and  eqn  36  can  now  be  re -expressed  as 

!(0+)  -  F(0" ) j  (39) 

Thus  we  see  that  the  change  in  the  initial  accelerations  is 
proportional  to  any  instantaneous  changes  (steps)  in  the 
magnitudes  of  the  applied  loads  at  t  =  0. 

To  determine  the  initial  accelerations  at  time  0  one  needs 
to  determine  the  accelerations  at  time  0  and  then  add  to  them  the 
additional  component  of  acceleration  from  the  change  in  load 
histories.  Should  there  be  a  continuous  transition  from  one  load 
history  to  the  other,  then 

u(0+)  =  u(0" )  (40) 

and  the  accelerations  at  time  0  are  the  accelerations  used  in  the 
initial  value  problem.  Satisfying  the  initial  conditions  on 
acceleration  in  this  manner  effectively  removes  the  a  order 
singular  terms  on  both  sides  of  eqns  32  and  33. 

The  remaining  singular  terms  in  these  equations  do  not  have 
corresponding  terms  on  the  force  side  of  the  equation.  To 
preserve  the  equality  one  must  conclude  that  the  coefficients  of 
these  singular  terms  are  zero.  Note  that  setting  these 

coefficients  to  zero,  in  effect  generates  the  remaining  initial 


u (0+)  -  u(O')  =  M' 


(37) 

(38) 
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conditions  needed  in  eqn  31.  From  eqn  32 


(m-2n-i)/3 


D 


u(0  ) 


£  =  1,2,3,  ■  , m- 2n- 1 


(41) 


and  from  eqn  33 


(m- 2n- i ) 


D 


u(0+) 


-  0 


£  -  1,2,3, •• • , m-2n- 1 . 


(42) 


Proof  is  given  in  Appendix  C.  Hence,  one  can  see  that  the 
initial  values  of  the  fractional  derivatives  of  displacement 
greater  than  second  order  and  less  than  order  /3m  must  be  zero  to 
preserve  the  equation  of  motion. 

Adding  the  two  equations  of  motion  and  recalling  that 


u(t)  +  u(t)  =  u(t) 


t  >  0 


(43) 


yields 

M ( 1  t  bDQ)u(t)  +  (k  +k  Da)u(t)  =  (l+bDQ)  ?(t)  +  G(t)  (44) 

=  =o  =1 

which  is  identical  to  eqn  8  except  for  one  very  important  detail: 
the  fractional  derivative  operator  has  changed  from  the  original 
definition,  eqn  1,  to  the  modified  definition  eqn  24.  Recall  that 
the  modified  basis  functions  use  this  modified  definition  as  well. 
In  fact  the  entire  initial  value  problem  (constituted  by  eqns  44, 
12,  25,  and  31)  and  its  solutions  (eqn  30)  can  be  cast  in  terms  of 
the  modified  definition  of  fractional  differentiation.  Using  the 
composition  property,  for  the  modified  operator, 


DQ[D7[u(t) j ]=DQ+7[u(t)] , 


(45) 
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which  only  applies  when  the  initial  values  of  the  fractional 
derivatives  of  the  displacements  are  zero,  one  can 
straightforwardly  demonstrate  that  eqn  44  leads  to  a  corresponding 
form  of  eqn  12  where  the  operator  is  replaced  by  .  Similarly, 

q  A  n 

the  Dpoperators  in  eqn  31  can  be  replaced  by  Dp  as  well.  Also 
noting  that  the  particular  solution  in  the  initial  value  problem 
is  in  effect  an  excitation  from  a  quiescent  state,  one  can 
demonstrate  that  the  first  m-1  terms  in  the  series  expansion  of 
the  kernel  in  eqn  30  add  out  when  the  modal  solutions  are  combined 
to  construct  the  particular  solution.  The  resulting  expression 
for  that  part  of  the  modal  solutions  which  describe  the  system 
response  is  given  below. 

y  (t)  =  y  (0)  E  f-(a.t)^ 

J  J  PK  J 

(46) 

E4'<v)1]vc‘,>dr 

Proof  is  given  in  Appendix  D.  The  modified  equality  symbol 
indicates  that  this  expression  is  true  only  in  the  context  of  the 
total  system  response. 

At  this  point  one  might  be  tempted  to  assert  that  the 
original  definition  of  fractional  order  differentiation,  eqn  1,  is 
somehow  inappropriate  for  the  initial  value  problem.  Not  true. 
Recall  that  the  initial  value  problem  has  insufficient  numbers  of 
physically  motivated  initial  values  to  determine  uniquely  the 
overall  homogeneous  solution  as  a  superposition  of  solutions  to 
the  modified  basis  equations.  The  additional  auxiliary  initial 
conditions,  developed  in  Appendix  C  by  suppressing  singular 


(-aV-1  D- 
1 ' 

o  J 
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behavior  at  time  zero,  provided  precisely  the  number  of  needed 


initial  conditions  for  a  unique  solution.  In  fact  the  original 
definition  generated  this  singular  behavior  without  which  the 
initial  value  problem  would  flounder  for  lack  of  initial 
information . 

To  test  the  robustness  of  this  generalized  initial  value 

problem,  one  needs  to  ascertain  its  ability  to  generate  the 

structural  response  to  impulsive  loading.  The  method  entails 

solving  the  initial  value  problem  for  a  step  response  (using 

initial  accelerations,  eqn  39)  from  a  quiescent  state  and  noting 

that  the  impulse  response  is  the  first  derivative  of  the  step 

response.  The  structural  response  for  a  unit  impulse  load  at  the 
til 

z  degree  of  freedom  of  the  structure  is 
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(97) 


where  z  is  an  N  order  column  vector  of  zeroes  except  the  z 
element,  which  is  one.  Again  the  solution  is  seen  to  be 
continuous  and  is  expressed  in  terms  of  the  modified  operator  and 
the  Mittag-Lef fler  function.  Derivation  of  this  expression  is 
given  in  Appendix  E. 

There  are  several  similarities  between  the  initial  value 
problem  posed  here  and  the  classical  initial  value  problem.  First 
the  equations  of  motion  are  seen  to  have  both  homogeneous  and 
particular  solutions  which  are  uniquely  dependent  on  the  initial 
values  and  the  forcing  functions  respectively.  In  addition 
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initial  values  for  all  rational  order  derivatives,  up  to  but  not 
including  the  highest  order  in  the  equation  of  motion,  are  needed 
to  determine  a  unique  solution.  Also,  the  impulse  response  is 
seen  to  be  a  homogeneous  solution  with  special  initial  conditions. 
These  similarities  suggest  that  other  parallels  do  exist. 
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Broader  Issues  in  the  Generalized  Initial  Value  Problem 

The  successful  construction  of  the  solution  to  the  rational 
order  initial  value  problem  yields  some  insights  into  the  nature 
of  the  underlying  mathematics.  One  of  the  more  unusual  results  is 

number  of  initial  conditions  needed  to  construct  a  unique 
solution.  It  appears  that  a  m/n  order  equation  needs  m  initial 
conditions:  one  condition  for  the  function  and  a  condition  for 
each  of  the  p/n  order  derivatives  p  -  1 , 2 , 3  ,  •  •  ■ m- 1 .  As  shown  in 
Appendix  C  the  physically  motivated  initial  conditions  -  initial 
displacement,  velocity  and  acceleration  -  are  complemented  by 
precisely  the  needed  number  of  auxiliary  initial  conditions. 
These  auxiliary  initial  conditions  are  the  result  of  strictly 
mathematical  considerations. 

Moreover,  how  does  posing  the  m/n  order  equation  as  a  2m/2n 
order  equation  change  the  nature  of  the  solution  for  which  2m 
initial  conditions  now  appear  to  be  needed?  In  fact  the  solutions 
in  both  cases  are  identical  from  which  one  infers  that  the  basis 

4 

functions  of  the  two  solutions  are  inter-related.  This  means  that 
the  1/m  order  and  l/2m  order  Mittag- Leff ler  functions,  eqns  28  and 
30,  are  related.  These  relationships  among  the  rational  order 
Mittag-Leff ler  functions  will  be  shown  later.  Choosing  to  solve 
the  equation  as  a  2m/2n  system  versus  a  m/n  system  leads  to 
higher  order  matrix  equations  and  simultaneously  provides  m  mere 
auxiliary  conditions  ensuring  a  unique  solution.  The 


One  needs  to  establish  that  different  but  apparently 
equivalent  forms  of  the  equations  of  motion  produce  identical 
solutions.  This  is  similar  to  posing  N  second  order  equations  as 
2N  first  order  equations  where  both  sets  of  equations  yield  the 
same  solution. 
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requirements  for  the  physically  motivated  initial  conditions 
remain  unchanged. 

Assume  that  the  largest  common  factor  for  the  fractional 
orders  is  1/n  where  the  modified  basis  equations  take  the  form 

(D1/n  +  aj/n)  y  j  ( t )  =  gj(t),  j  -  1,2,3,  •  •  •  ,  m  (48) 

but  one  chooses  to  solve  the  problem  with  a  basis  fraction 
of  l/2n . 


/rs2 /2n  l/2n.  ,  .  ,  x 

(D  +  »p  )y  (t)  -  gp(t) 


p  -  1,2,3, -  ,2m  (49) 


For  each  basis  equation  in  eqn  48  with 

there  are  two  corresponding  equations 

.  /  1/n  ,  .  l/2n  , 

±v-a.  ,  or  ±ia.  .  The  solutions 

J  J 

of  eqn  49,  in  the  form  of  eqn  30,  are 


characteristic  value  a. 

J 

in  eqn  49  with  roots  of 
to  these  two  equations 


yp(t)  -  yp(0)  E1/2n 


■iU.t,1'2"] 
J  > 


D1 - l/2n 

E 

[  -  i(a  .  r )  1/2n]"| 

l/2n 

l  J  JJ 

g  (t-r)dr 
P 


(50) 


and 


Vl(t)  =  yp+l(0)  El/2n 


Vi(a.t)1/2^ 


yC 


D 


1  -  l/2n 


'l/2n 


l/2n\  -i 


+  i ( aj  f  ) 


Jgp+1(t-r)dr 


(51) 


Solving  the  problem  with  a  basis  fraction  of  l/2n  instead  of  1/n 
produces  m  additional  auxiliary  initial  conditions  of  the  form 
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Dp/2n  u(0) 


p  odd  and  0  <  p  <  2m 


(52) 


The  solutions  to  eqn  48  also  satisfy  these  initial  conditions.  In 
Appendix  C  it  is  shown  that  if  u(0)  exists,  then  the  initial  value 
of  any  rational  order  derivative  less  than  one  is  zero,  while  u(0) 
may  be  independently  specified. 

The  solutions  given  in  eqns  50  and  51  for  eqn  49  when 
combined  should  produce  a  corresponding  solution  to  eqn  48,  which 
again  based  on  eqn  30,  takes  the  form 


y.(t)  -  yj(0)E1/n 


’-U.t,1/”] 

l  J 


Dl-1^ 

E  . 

[  -  (a  .  r  ) 1/n) 1 

.  1/n 

13  JJ 

g  ( t  -  r ) dr 


(53) 


However,  the  solutions  in  eqns  50  and  51  are  not  combined 
arbitrarily.  The  solutions  are  combined  in  a  manner  consistent 
with  solving  eqn  48  as  a  system  of  twice  as  many  equations  of  half 
the  order  of  differentiation.  As  such,  the  loads  g  (t)  may  be 
viewed  as  modal  loads  generated  bv  g.(t),  and  applying  the 
following  conditions 

y. (0)  =  yp(0)  +  yp+1(0)  (54) 

D1/2V(0)  -  0  (55) 

the  resulting  expression  for  v.(t)  is 
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where  the  integral  terms  are  expressed  in  terms  of  eqn  30. 
However,  defining  the  functions  in  the  brackets  as  generalized 
cosine  and  sine  functions 


Dl/2nsini/2n(ajt)  -  a.1/2ncos1/2n(a.t)  (60) 

Dl/2ncoSi/2n(ajt)  -  -a.1/2nsin1/2n(a.t)  (61) 
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More  general  properties  of  the  fractional  order  sine  and  cosine 
functions  and  their  hyperbolic  counterparts  are  given  in  Appendix 


F.  The  appearance  of  generalized  sine  and  cosine  functions  can  be 
explained  by  observing  that  eqn  48  could  have  been  equivalently 
posed  as  a  generalized  second  order  differential  equation  of  basis 
l/2n . 


(D2/2"  + 


2/2n 


)yj(t)  - 


gj  (t) 


(62) 


To  complete  the  development  of  the  relationship  between  the 
basis  1/n  and  l/2n  functions,  the  composition  property,  eqn  7,  is 
applied  to  the  derivative  operator  in  the  integrand  of  eqn  59 


D1  -  X/2n 


1  -  1/n  +  l/2n 


(63) 


and  noting  that 


nl/2n  . 

D  SLnl/2n(ajt) 


l/2n 

aj  C°Sl/2n(ajt) 


(64) 


eqn  59  now  becomes 


y.(t)  -  y.(0)cos1/2n(a.t)  + 


rc 

D 


1  -  1/n 


cosl/2n(ajr) 


g(t-r)dr  (65) 


Comparing  this  solution  to  the  solution  of  eqn  48  given  in  eqn 
53,  one  concludes  that 


El/nl  (ajt} 


1/nl 


cos.,  (a.t) 
l/n  J 


Jl/2n 


i(a1t)1/2n)'El/2n 


■  i (3j  t) 


l/2n 


(66) 


or 


2E 


l/n 


l/n 


l/2n) 


+  E 


l/2n 


- i(a^  t) 


l/2n'l 


(67) 


For  aj’/'n  negative  the  relationship  becomes 


[(a.t)1/"] 

=  E 

f(a.t)1/2"1 

+  E- 

[-(a.t)1/2"] 

l  J  J 

l/2n 

[  J 

l/2n 

l  J 

"  2  COShl/2n(ajt) 


(68) 


These  relationships  can  be  straightforwardly  verified  by  expanding 
the  Mittag-Leff ler  functions  in  terms  of  their  series  definition, 
eqn  26,  and  adding  terms  of  equal  powers.  The  more  general 
relationship,  also  verified  in  the  same  manner,  takes  the  form 


"m/n 


±(at) 


m/n 


1 

P 


m/pn 


(69) 


where  (-l)^^  are  the  p  different  1/pth  roots  of  minus  one.  This 
relationship  allows  one  to  demonstrate  the  equivalence  of 
solutions,  u(t) ,  when  determined  in  terms  of  different  basis 
fractions . 

One  now  concludes  that  regardless  of  the  choice  of 
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appropriate  basis  fraction,  /3 ,  and  the  representation  of  the 
solution,  the  uniqueness  of  the  solution  may  be  demonstrated  using 
eqn  69.  One  also  concludes  that  a  change  in  the  basis  fraction 
adjusts  the  number  of  auxiliary  initial  conditions  to  preserve  the 
uniqueness  of  the  solution,  leaving  the  role  of  physically 
motivated  initial  conditions  unaltered.  Moreover,  the 

relationships  among  different  order  basis  function  leads  to 
generalized  definitions  of  sine,  cosine  and  the  hyperbolic 
functions  (Appendix  F)  .  The  observation  that  these  generalized 
functions  degenerate  into  the  definitions  for  the  normal  sine, 
cosine,  etc.,  when  fractions  are  set  to  one,  coupled  with  the 
observation  that  eqn  53  becomes  the  general  solution  to  a  first 
order  equation,  eqn  62,  when  n  is  set  to  one,  lends  credence  to 
the  view  that  these  fractional  order  differential  equations  and 
their  solutions  are  legitimate  generalizations  of  their  integer 
order  counterparts.  Furthermore,  this  solution  technique  may  be 
viewed  as  a  legitimate  generalization  of  the  initial  value 
problem . 

Yet  another  similarity  with  ordinary  differential  equations 
arises  when  one  examines  the  nature  of  homogeneous  solutions  for  a 
fractional  order  differential  equation  having  repeated  roots  in 
its  characteristic  equation.  Posing  the  differential  equation  as 
a  system  of  /3  order  modified  basis  equations  in  Jordan  form 
(10:35),  the  two  basis  equations  containing  the  same  root  take 
the  form 
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(70) 


D^y.  +  sP.  y.  -  0 

J  J  J 


D^y .  +a^y.  ..  +  y .  =0 

^J  +  l  J  7J+1  yJ 


(71) 


Solving  for  y^  (eqn  28)  and  using  y.  as  a  forcing  fur.ction  to 


determine 


y^+^  (eqn  30)  yields 


^j+l(t)  “  yj+l(0)  E/3 


-(a  t/1  +  y .  (0)  /3_1t  D1-/? 

1 

m 

3 

l  J  J  J 

L  P\ 

l  J  JJ 

(72) 


V  .en  /9  is  one  this  solution  becomes  the  solution  to  a  second  order 
differential  equation  with  repeated  characteristic  values. 


y-j+x^)  -  yj+i<°>  e'aJt  +  t  yj(0)  e'ajt 


(73) 


When  the  root  is  repeated  more  than  once,  y  ^(t)  is  used  to  find 
yj+2(c)  and  so  on  until  the  solution  is  complete. 


+  ,0.  yJ+2(t)  +  y.+1  -  0 


(7d) 


The  general  form  of  the  homogeneous  solution  for  m  repeated  roots 


is 


y  (t) 
j+m- 1 


V0)  E/? 


~(ajc) 


m_  1  rk-i  r  i 

X  yi+k  i(°)  n  (/9p)’tDi^ 

k=2  J  Lp=1  l  J 


“(a  t) 


P) 


(75) 


which  for  /3  =  1  becomes  the  solution  for  an  ordinary  differential 
equation  having  constant  coefficients  with  m  repeated  roots  in  its 
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characteristic  equation. 

When  one  is  confronted  with  all  of  the  similarities  between 
fractional  order  and  ordinary  differential  equations,  it  seems 
appropriate  to  examine  the  fractional  order  initial  value  problem 
as  a  well  posed  problem.  The  existence  of  solutions  has  been 
established.  The  uniqueness  of  solutions  can  be  readily- 
established  by  limiting  the  class  of  loads  and  responses  to  those 
having  Laplace  transforms.  One  can  expand  this  class  of  functions 
to  include  those  loads  that  produce  unique  responses  through  the 
use  of  the  Green's  functions,  eqn  46  and  47,  and  convolution. 
Continuous  dependence  on  the  data  can  be  established  using  eqns 
28,  31,  46,  and  47  and  noting  that  the  Green's  functions  are 
non-singular.  In  effect  the  fractional  order  differential 
equations  exhibit  precisely  those  characteristics  expected  from 
ordinary  differential  equations. 

One  might  ask  why  not  use  Laplace  transforms  to  solve  the 
problem  and  dispense  with  the  formalism  constructed  here.  First, 
one  is  not  now  restricted  to  those  loading  and  displacement 
histories  having  Laplace  transforms,  although  in  practice  most 
loadings  and  responses  of  engineering  interest  do  have  Laplace 
transforms.  In  principle  it  is  possible  to  formulate  the  initial 
value  problem  with  Laplace  transforms  using  delays  to  stop  and 
start  loading  histories,  but  several  key  features  of  the  initial 
value  problem  are  missed.  The  transition  to  the  modified  operator 
and  the  rationale  for  the  formulation  of  the  auxiliary  initial 
conditions  can  be  completely  overlooked.  Even  if  the  modified 
operator  were  discovered  using  Laplace  transforms,  the  initial 
values  of  the  fractional  derivatives  do  not  appear  in  the 
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formulation,  because  transform  theory  has  correctlv  set  them  to 
zero.  The  fact  that  they  do  not  appear  obscures  their  existence 
and  the  existence  of  the  auxiliary  conditions  needed  for  a  unique 
solution.  As  a  direct  consequence,  the  non-singular  nature  of  the 
Green's  function  solutions  is  missed  as  well.  Also  missed  are  the 
generalized  nature  of  the  basis  functions  being  special 
Mi ttag- Lef f ier  functions,  the  inter-relationships  among  the 
different  sets  of  basis  functions  and  the  definitions  of  the 
generalized  sine,  cosine,  etc.,  functions  spawned  by  these 
inter-relationships  . 

In  essence,  Laplace  transforms  by  themselves  do  not  readily 
reveal  the  structure  of  the  initial  value  problem,  and  the  needed 
auxiliary  initial  conditions  for  a  unique  solution  are  not 
produced,  except  possibly  in  retrospect.  As  a  result  the  major 
features  of  this  generalized  initial  value  problem  can  be  missed, 
and  one  is  never  justified  in  claiming  that  the  class  of 
differ- integral  equations  solved  here  are  in  fact  generalized 
diffe^entiai  equations. 
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Conclusions 


The  most  encompassing  conclusion  is  that  the  class  of 
differ- integral  equations,  represented  by  eqn  8,  may  be  viewed  as 
a  generalization  of  ordinary  differential  equations  with  constant 
coefficients.  When  an  equation  of  motion  is  posed  in  terms  of  m 
differential  equations  of  fractional  order  j8 ,  m  sets  of  initial 
conditions  are  needed  to  determine  a  unique  solution.  Careful 
examination  of  the  singular  behavior  generated  by  fractional 
differentiation  leads  to  auxiliary  initial  conditions  that 
supplement  the  physically  motivated  initial  conditions  producing  a 
total  of  m  sets  of  initial  conditions  and  a  unique  solution.  The 
resulting  fractional  order  differential  equations  are  posed  in 
terms  of  a  modified  definition  of  fractional  differentiation, 
which  appears  to  act  on  the  solutions  as  though  they  were 
analytically  continued  from  negative  time  into  positive  time. 
This  is  precisely  the  mathematical  property  needed  for  a 
generalized  initial  value  problem. 

Several  other  key  features  of  the  mathematical  development 
ler.J  credence  to  the  claim  that  the  mathematics  herein  embodies  a 
generalization  of  ordinary  differential  equations.  In  every 
instance,  setting  the  basis  fraction  to  one  transforms  the 
solutions  into  the  traditional  solutions  associated  with  ordinary 
differential  equations.  In  addition  the  solutions  to  the  modified 
basis  equations  are  posed  in  terms  of  Mittag-Leffler  functions, 
long  viewed  as  the  fractional  order  generalization  of  the 
exponential  function.  The  general  form  of  the  basis  solutions 
contains  a  portion  dependent  on  the  initial  condition  and  a 
portion  dependent  the  forcing  function,  clearly  identifiable  as  a 
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generalized  homogeneous  solution  and  a  generalized  particular 
solution,  respectively.  The  form  of  the  homogeneous  basis 
solution  for  the  case  of  repeated  character stic  roots  is  strongly 
reminiscent  of  its  integer  order  counterpart,  and  the  two 
expressions  are  identical  when  the  basis  fraction  is  one.  The 
generalized  sine  and  cosine  functions  along  with  the  generalized 
hyperbolic  functions  are  also  very  similar  to  the  traditional 
functions.  These  similarities  lead  to  generalized  identities  that 
also  become  traditional  trigonometric  identities  when  the  basis 


fraction  is  one.  In 

addition  to 

all 

these  similarities 

with 

ordinary  differential 

equations , 

the 

observation  that 

the 

fractional  order  initial  value  problem  appears  to  be  a  well-posed 
problem  makes  the  case. 

This  evidence  leads  one  to  conclude  that  the  strong 
similarities  between  the  generalized  initial  value  problem  and  the 


traditional 

initial  value 

problem 

are  not 

coincidence,  and  that 

this  generalized  initial 

value 

problem 

may 

be  viewed  as  a 

legitimate 

extension  of 

the 

theory 

for 

ordinary  linear 

differential  equations  with  constant  coefficients. 
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APPENDIX  A  -  Differentiating  When  Time  is  Negative 

This  section  presents  the  technique  needed  to  perform 
fractional  differentiation  over  the  time  interval  prior  to  the 
initial  time  of  the  initial  value  problem.  In  the  process  of 
setting  up  the  initial  value  problem  one  needs  to  examine  the 
mathematical  behavior  of  the  system  just  prior  to  and  just  after 
the  initial  time  (t  =  0)  to  determine  the  nature  of  additional 
initial  conditions. 

Unfortunately,  the  definition  of  fractional  order 
differentiation,  eqn  1,  becomes  ambiguous  when  the  lower  limit  of 
integration  is  taken  as  a  negative  value. 


dq  [y(t)] 


r(i-a)  dt  J 


y(u) 

(t-u)Q 


du 


(A  - 1 ) 


For  t  less  than  zero  and  greater  than  -tQ,  the  kernel  in  the 
integral  becomes  a  complex,  multivalued  relationship.  The 
definition  of  fractional  differentiation  given  in  eqn  1  does  not 
readily  allow  one  to  specify  time  zero  strictly  as  a  matter  of 
convenience.  It  appears  that  this  definition  of  fractional 
differentiation  insists  that  time  zero  be  chosen  at  times  prior  to 
any  non-zero  values  of  y(t)  for  the  definition  to  produce  at 
unique  result. 

One  would  like  to  have  the  latitude  of  choosing  time  zero  to 
coincide  with  the  initial  time  of  the  initial  value  problem.  To 
circumvent  the  ambiguity  raised  in  eqn  A-l,  one  may  employ  the  two 
time  scales  shown  in  figure  1  on  page  15:  the  t  scale  to  set  up 
the  overall  problem  and  the  t  scale  to  accommodate  the  initial 


value  problem.  The  first  time  scale,  t,  takes  zero  to  be  the 


the  interval  t  )  in  terms  of  t  produces  the  equation  for  prior 
motion 


M  (l+bDQ)  u  (t  -  t  )  +  (k  +  D“)  u  (t  -  tQ) 


Q:  v  —  - 


(A-5) 


=  ( 1+bD 


a)  jf(t  - 


t  )  -  U(t  -  t  )F(t  -  t  ) 

O  0—0 


Using  eqn  A~4  to  construct  these  derivatives,  then  applying 
Leibnitz's  rule  for  differentiating  the  integrals  to  separate  out 
the  singular  behavior  and  substituting  t  for  t  —  t  produces  eqn 
32.  Those  terms  resulting  from  this  process  not  specifically 
identified  in  eqn  32  constitute  G(t)  ,  the  pseudo  forcing 
functions,  which  produces  the  residual  motion,  u,(t  >  0)  due  to 
previous  loading  (t  <  0)  of  the  structure.  This  residual  motion 
is  superimposed  on  the  response  of  the  structure,  u(t) ,  to  present 
loading,  F(t),  to  produce  the  total  response  of  the  structure  for 
t  >  0. 
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onset  of  any  recent  history  of  motion  for  the  structure  and  the 
second  scale,  t,  takes  zero  to  be  the  later  time  about  which  one 


wants  to  delineate  previous  and  present  loading  histories.  This 

later  time,  t,  is  the  time  scale  of  the  initial  value  problem 

where  the  initial  time  is  t  =  0.  The  only  difference  in  the  time 

scales  is  the  shift  factor  t  ,  the  time  interval  of  the  recent 

o 

history  of  motion.  Hence  the  relationship  between  the  two  time 
scales  is 

t  -  t  —  t  ( A  -  2 ) 

o 

Posing  eqn  A-l  in  terms  of  this  time  shift  yields 


y(t-to) 


1  d 


=r 


y(u) 


r(1->  dtJ-t  (t  -  t  -  u)“  du 

0  O 


(A-  3) 


and  with  the  change  of  variable  u  -  T_t0  this  relationshop  becomes 


y(t-to) 


m-«) 


y 


t  y(r-t  ) 


dr 


dt  J  0  (t-r) 


(A-4) 


which  appears  in  eqn  17,  with  a  set  equal  to  /S . 

The  shift  of  t  in  the  time  scale  restores  the  sinzle  valued 
o  ° 

property  of  the  fractional  derivatives  operator.^  More 
importantly  eqn  A-4  allows  one  to  evaluate  the  behavior  of 
fractional  derivatives  in  the  equation  of  motion,  eqn  7,  prior  to 
the  initial  time  of  the  initial  value  problem  which  occurs 
for  t  =  t  Posing  the  equation  of  motion  for  prior  loading  (over 


Note  that  replacing  t  with  t  and  setting  t  =  0  produces  the 
original  definition,  eqn  1. 


APPENDIX  B  -  The  Structure  of  the  Eigenvector 

This  section  establishes  the  systematic  structure  of  the 
eigenvectors  (the  columns  of  the  orthonormal  transformation  matrix 
in  eqn  14)  associated  with  the  expanded  equations  of  motion,  eqn 
12.  In  particular  one  demonstrates  that  the  j  —  eigenvector  eqn 
12  has  the  following  structure 
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where  is  the  corresponding  j —  eigenvector  associated  with  the 
homogeneous  form  of  eqn  11  and  -sP  is  the  eigenvalue  associated 
with  both  eqns  11  and  12.  Staged  more  succinctly,  one  will  prove 
that 
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given  the  homogeneous  form  of  eqn  12  shown  below 
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To 

begin  the 

proof  one  focuses  on  the 

upper 

m-1  sets  of 

matrix 

equations 
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re-expressed  as 
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The  absolute  value  of  the  determinant  of  [A],  indicated  by  det[A] 
can  be  expressed  as 


m- 1 


det[A]  =  b  (det  M) 


,  m- 1 


(B-6) 


where  M  is  the  system's  mass  matrix.  Since  the  mass  matrix  is 
positive  definite,  then  the  determinant  of  [A]  is  always  non— zero 
and  [A]  is  not  singular.  Hence  the  only  solution  to  eqn  B-5  is 
the  trivial  solution,  and  as  a  result 
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*3j  *  V-?> 
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from  which  one  concludes 
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APPENDIX  C  -  Deriving  the  Auxiliary  Initial  Conditions 

This  section  establishes  the  sources  of  the  auxiliary  initial 
conditions.  These  auxiliary  conditions  are  needed  to  determine  a 
unique  solution  for  the  initial  value  problem  when  the  equations 
of  motion  are  posed  as  a  system  of  /5  order  differential  equations. 

The  auxiliary  conditions  have  two  sources.  The  first  is  eqn 
33,  the  equation  of  motion,  where  certain  singular  terms  appear  on 
the  response  side  of  the  equation  without  corresponding  singular 
terms  appearing  in  the  applied  forces.  The  equation  of  motion  is 
preserved  by  setting  the  coefficients  of  these  singular  terms  to 
zero.  These  coefficients  are  the  initial  values  of  the  fractional 
derivatives  of  order  greater  than  2  and  less  that  2  +  a  of  the 
structural  displacement  histories.  The  source  of  the  remaining 
auxilary  conditions ,  those  derivatives  of  order  between  one  and 
two  and  between  zero  and  one,  is  the  following  theorem.  Given  a 
set  of  functions,  u(t),  which  are  a  linear  combination  of  Mittag- 
Leffler  functions  as  in  eqn  31  where  u(t)  have  bounded  first 
derivatives  at  t  =  0 ,  then  all  fractional  derivatives  of  these 
functions  of  rational  order  less  than  one  have  initial  values  of 
zero  at  t  =  0. 

The  proof  starts  with  eqn  32. 
m  •  N 

^(t)  =  l  ix.  y.(t)  (32) 

j=l  J  J 


where  <f>  are  the  eigenvectors  associated  with  eqn  11  and  y.(t) 

— J 

are  Mittag-Lef f ler  homogeneous  solutions  for  the  modified  basis 
equations  shown  in  eqn  28.  Identifying  these  homogeneous  solutions 
appearing  in  eqn  31  with  the  subscript  j  produces 


yj(t)  -  yj(0)  Ef 


k,  J 


(C-1) 


and  using  eqn  27  Co  express  these  functions  in  series  form 
produces 


■ (  a  j  t ) 


-  X 


f— (a.t)'l 
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J 
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P=o  r(i+p/3) 

Substituting  eqn  C-2  into  eqn  C-1  and  then  substituting  the 
resulting  equation  for  y  (t)  in  eqn  32  yields 


m  •  N 


Vc)  “  .^ij  yj(0)  J0  mi^) 
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Taking  the  first  time  derivative  of  this  expression  results  in 


m  ■  N  ® 

u(t>  =  X  h*  y/°)  X 

j=l  J  J  P-1 


(~apP  p/3  t 
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where  the  derivatives  of  the  constant  terms  in  the  infinite  series 
(p=0)  are  now  zero.  Because  one  of  the  two  summations  has  a 
finite  number  of  terms,  the  order  of  summation  may  be  interchanged 
without  loss  of  generality.  Also  noting  that 


P  £  1 

r(l+p/3)  "  r(p/3) 

the  expression  for  the  derviative  u,  (t)  becomes 


Vc) 


«  p/3-1 

P?i 
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y.  (0) 
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Recall  that  u(0)  is  assumed  bounded  and  at  worst  has  a  step 


discontinuity  at  t  -  0,  and  that  /3  is  of  the  form  1/n.  Therefore, 
the  first  n-1  terms  of  the  infinite  series  must  be  zero  because 


they  are  singular  and  unbounded  at  t=0 .  Setting  the  coefficients 
of  the  first  n-1  terms  equal  zero  yields 


m-N  R  p 

Y  4.  .(-aP)  y.(0)  -  0  p  =  1,2,3, ■ ,n-l 

j L1  -lj  J  J 


Applying  the  results  of  Appendix  B 


k  =  1-2,3,  •••  ,m 

to  equation  C-7,  while  noting  the  m  >  2n,  results  in 


( C  -  7  ) 
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1  *  ,  7,(0)  -  0  p  -  1,2,3, ••• ,n-l  (C-8) 
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Examining  the  structure  of  the  j —  eigenvector  for  eqns  12, 
shown  in  equations  B-l,  and  using  eqn  31  one  concludes  that 

m  ■  N 

°P/r‘  V0)  =  £  £pl  yi(0)  P  -  1,2,3,  •••  ,n-l  (C  -  9 ) 

j=l  J 

and  with  eqn  C-8,  one  can  further  conclude  that 


DP/n  Uh(°)  =  0  p  -  1,2,3, ,n-l  (C-10) 

Because  n  has  not  been  specified,  in  principle  any  n  could 
bechosen.  One  can  conclude  that  p/n  can  take  on  any  rational 
value  between  zero  and  one  and  the  proof  of  the  theorem  is 
complete . 
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However,  the  theorem  need  now  be  applied  to  determine  the 

initial  values  of  the  derivatives  of  u^(t)]^  of  rational  order 

between  one  and  two.  The  most  intense  forces  applied  to  the 

structure  in  the  initial  value  problem  are  the  step  functions 

turning  off  and  turning  on  the  loading  histories  at  time  zero.  A 

step  load  is  incapable  of  instantaneously  changing  the  system's 

kinetic  energy  and  one  concludes  that  the  resulting  velocity 

histories  are  continuous  functions  for  all  time.  Step  loading  is 

also  incapable  of  producing  unbounded  accelerations,  therefore  the 

first  derivative  of  velocity  is  bounded  and  piecewise  continuous. 

From  this,  one  can  conclude  that  the  velocities  v,  (t)  have  bounded 

— h 

derivatives,  v  (t)  ,  and  the  above  theorem  applies  to  the  initial 
values  for  the  fractional  derivatives  of  rational  order  between 
zero  and  one  of  the  velocities.  The  result  is 

Dp/n  vh(0)  =  0  p  -  1,2,3, ••• ,n-l  (C-ll) 

Expressed  in  term  of  the  displacements,  these  conditions  are 

Dp/n+1uh(0)  =  0  p  =  1,2,3,  •••  ,n-l  (C-12) 

and  are  seen  to  be  the  initial  values  of  the  derivatives  of 
rational  order  between  one  and  two  of  the  displacements. 

Given  that  the  initial  displacements  and  the  initial 
velocities  are  known,  coupled  with  the  auxiliary  conditions  given 
in  eqn  C-10  and  C-12,  the  only  remaining  initial  conditions  needed 
in  eqn  31  are  the  initial  accelerations  and  the  initial  values  for 
the  fractional  derivatives  of  u^(t)  of  rational  order  between  2 
and  2+q.  These  conditions  are  determined  by  examining  the 
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singular  behavior  in  the  equation  of  motion  for  previous  loading, 
eqn  32,  and  the  equation  of  motion  for  present  loading,  eqn  33. 

The  singular  behavior  in  these  equations  of  motion  is 
generated  by  the  a  order  derivatives  of  the  step  functions  turning 
off  and  turning  on  the  previous  and  present  loading  histories  and 
by  the  2+a  order  derivatives  of  the  respective  resulting 
dispacement  histories.  The  singular  behavior  of  the  forcing 
function  and  the  a  order  derivative  of  the  displacements  in  the 
equation  of  motion  for  present  loading,  eqn  33,  is  derived  using 
eqn  23.  Calculating  the  2+a  derivative  of  the  displacements  is 
relatively  straightforward 


^(t)  -  ^(0)  +  u^(0)  •  t  +  X 

p=0 


2+p/S  m-N  p+2n 

£  V"?  yj<0)  <c'13) 


The  displacements,  shown  here  with  the  initial  displacements  and 
velocities  specified  and  the  auxiliary  conditions  of  eqn  C-10  and 
C-12  applied,  are  first  differentiated  twice  and  then  successive  / 3 
order  derivatives  are  taken  to  produce  the  2+a  order  derivative. 
The  result  is 


D2+Q  u^(t)  =  l 

p=0 


tP/3-a 
T ( l+p/3-a) 


m-N  p+2n 

I  y,(0) 


(C- 14) 


and  recalling  the  /?  =  1/n  and  from  eqn  10  that  a  -  /Jq ,  this 
expression  becomes 
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Separating  out  the  singular  terms  yields 
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For  p  -*  0  the  order  of  the  singular  term  is  a  and  the  coefficient 
of  this  singular  term 


m-N  2n 

ij(-aj}  yj(0)=Sh(0+) 


(C-17) 


appears  in  eqn  44  as  the  initial  value  of  the  acceleration  due  to 
the  present  loading. ^  Recall  that  eqn  35  resulted  from  setting 
equal  the  coefficients  of  the  a  order  singular  terms  on  both  sides 
of  eqn  33. 

The  remaining  singular  terms  in  eqn  C-16  do  not  have 
corresponding  singular  terms  in  the  forces.  To  preserve  the 
equality  of  eqn  33  one  must  conclude  that  the  coefficients  of 
these  singular  terms  are  zero.  Setting  these  coefficients  to  zero 
yields 

m-N  p+2n 

I  i ( )  y,(0)-0  p  -  1,2,3,-. -.q-l  (C - 18 ) 

j-1  J  J  J 

Again,  using  the  results  of  Appendix  B,  eqn  B-2,  eqn  C-18  becomes 


This  result  can  be  confirmed  by  differentiating  C-13  twice  and 
setting  t  -  0. 
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m  •  N 

l  t  y.(0+)=0  £  =  2n+l,2n+2,2n+3, • • • ,m-l  (C-19) 

j-1  J  J 

This  equation  sets  the  initial  values  of  the  fractional 
derivatives  of  the  displacements  of  rational  order  between  2  and 
2  +  q  to  zero. 

Applying  the  process  that  led  to  eqns  C-14  through  C-19  to 
the  equation  of  motion  for  previous  loading,  eqn  32,  and  equating 
the  a  order  singular  terms  leads  to  eqn  34.  The  methods  used  to 
determine  the  fractional  derivatives  at  time  0  are  presented  in 
Appendix  A.  Combining  eqns  34  and  35  produces  the  initial 

condition  on  acceleration,  eqn  39,  and  with  the  initial  values  of 

✓ 

displacements  and  velocities  completes  the  needed  physically 
motivated  initial  conditions.  These  conditions  coupled  with  the 
auxiliary  initial  conditions  established  in  eqns  C-10,  C-12  and 

C-19  constitute  the  complete  set  of  initial  conditions  needed  to 
determine  a  unique  solution  for  the  initial  value  problem. 
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APPENDIX  D  -  The  Basis  Green ' s  Function 

The  overall  particular  solution  (forced  response)  for  the 
equations  of  motion,  eqn  8,  can  be  constructed  from  the  solutions 
of  the  modified  basis  equations  appearing  in  eqn  30.  The  overall 
solution  is  a  linear  combination  of  the  modified  basis  particular 
solutions  prescribed  by  the  orthonormal  transformation  for  the 
equations  of  motion.  These  linear  relationships  are,  of  course, 
the  analogue  of  eqn  31  applied  to  the  particular  solutions. 
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Although  one  need  only  employ  the  lowest  set  of  matrix  equations 
to  determine  the  particular  response,  this  relationship  is  useful 
in  demonstrating  that  the  overall  particular  solution  does  not 
contain  the  response  generated  by  the  singular  terms  in  the  kernel 
appearing  in  eqn  30.  The  kernel  is  of  the  form 
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These  kernels  are  the  particular  solutions  when  the  forcing 

functions,  g  (t),  are  unit  impulse  functions. 

To  examine  the  behavior  of  the  impulse  response ,  one 

calculates  first  the  response  to  a  unit  step  load  and  then  takes 

the  first  derivative  of  the  step  response  to  produce  the  impulse 

cVi 

response.  The  response  of  the  j~  modified  basis  equation  for  a 
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unit  step  applied  at  the  z —  degree  of  freedom  is 
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Here  is  the  transpose  of  the  j~  eigenvector  of  eqn  11,  z  is  a 


th 


column  matrix  of  zero  elements  except  the  z —  which  is  one  and 
U(t)  is  the  unit  step  function.  Substituting  eqn  D-4  into  eqn 
D-3  and  evaluating  the  integral  produces  the  step  response, 
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I  is  the  fractional  integral  of  order  /3  and  is  defined  as 
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The  /3  order  integral  of  the  /?  order  Mittag- Lef f ler  function 
appearing  in  equation  D-5  can  be  shown  to  be 
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Using  eqns  D-5  and  D-7  and  summing  all  the  step  responses  of 


the  basis  functions  using  the  lowest  set  of  matrix  equations  in 
D-l  produces  the  structure's  step  response. 
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However,  noting  that  all  the  basis  step  responses,  eqn  D-5,  are 
zero  at  t  =  0  (see  eqn  D-7),  one  can  conclude  from  eqn  D-l  that 
all  of  the  fractional  derivatives 


D^ku  (t)"|  -0  k  -  0,1, 2, 3, • • • ,m-l  (D-9) 

Sz  t-0 

of  the  step  response  are  zero  at  t  =  0  as  well.  In  other  words, 
if  all  the  step  responses  are  zero  at  t  «  0,  then  all  of  the 
elements  in  the  right  column  vector  of  eqn  D-l  are  zero  for  t  =  0 
and  eqn  D-9  follows. 

It  follows  from  applying  eqn  D-9  to  eqn  D-8  that  the  first 

m-1  terms  in  the  time  series  representation  of  u  (t)  ,  eqn  D-8, 

z 

must  also  be  zero.  The  resulting  expression  for  the  step  response 
is 
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Taking  the  first  derivative  of  u  produces  the  impulse  response, 

z 
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u  (t).  Recall  that  /3  -  l/n. 
—o 
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m-N  _  -  (-«62n+tl+i-1  t1+<>+^ 

-«z<t:i  '  ^lj  hj  i  1  F(2 (D-12) 

However,  Che  overall  impulse  response  may  also  be  expressed  as  a 
linear  combination  of  the  effective  impulse  responses  associated 
with  each  mode  as  shown  below. 

m  ■  N 
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Using  eqns  D-12  and  D-13  one  can  determine  the  expression  for  the 
effective  impulse  response  associated  with  each  mode,  and  it  is 
seen  to  be 
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Given  the  effective  modal  impulse  response  one  can  now  work 
back  to  determine  the  new  form  for  the  kernel  of  the  convolution 
integral  in  the  particular  solution.  The  new  kernel  is 
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The  first  term  in  the  kernel  is  order  1+a .  The  earlier  expression 
for  the  kernel,  eqn  D-2,  shows  a  first  term  of  order  /3— 1 .  The 
difference  in  these  expressions  is  based  on  the  observation  that 
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eqn  D-9  demonstrates  that  the  first  m-1  terms  of  the  basis  step 
responses  sum  to  zero  when  constructing  the  overall  step  response 
using  eqn  D-l.  Consequently,  the  first  m-1  terms  of  the  basis 
impulse  responses  sum  to  zero  when  constructing  the  overall 
impulse  response.  Without  loss  of  generality  these  m-1  terms  may 
be  left  out  of  the  kernel  when  calculating  the  overall  particular 
solution  for  any  loading  history.  The  general  form  of  the  overall 
particular  solution  is  given  in  eqn  46. 
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APPENDIX  E  -  The  Response  to  Impulsive  Loading 

In  Appendix  D  the  impulse  response  or  Green's  function  was 
based  on  g  (t),  eqn  22,  being  impulsive  where  g^(t)  is  a 
combination  of  the  pseudo  forces  and  the  stress  operator  from  the 
viscoelastic  model  acting  on  the  applied  forces.  This  appendix 
presents  the  general  form  of  the  response  of  the  structure  to  an 
impulsive  force.  The  derivation  is  based  on  the  step  response  for 
which  the  foundation  has  been  laid  in  eqns  32  through  44  and 
Appendices  A,  B,  and  C.  Paralleling  Appendix  D,  the  determination 
of  the  impulse  response  rests  on  the  observation  that  in  all 
linear  systems  the  first  derivative  of  the  step  response  is  the 
impulse  response . 

The  derivation  of  the  step  response  is  based  on  solving  eqn 
44.  Note  that  for  step  loading  with  no  previous  motion,  G(t)  and 

A 

D  F(t)  are  zero  in  eqn  44.  Casting  the  resulting  equation  of 
motion  in  expanded  format,  eqn  12,  allows  one  to  solve  for  the 
step  response  by  specifying  the  initial  accelerations  in  the  state 
vector.  The  expressions  for  the  initial  accelerations  based  on 
eqn  39  is 

u(0+)  =  M'1|(0+)  (E-l) 

where  the  tildes  will  be  dropped  because  there  is  no  longer  a  need 
to  distinguish  between  previous  and  present  loading  and  responses. 

The  remaining  force  vector  in  eqn  44,  F(t),  is  taken  to  be  a 
unit  step  load  applied  at  the  z~  degree  of  the  freedom.  This 
load  vector  z  is  all  zeros  except  for  the  z~  element,  which  is 
one.  The  resulting  expression  for  the  initial  accelerations  is 
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u(0+)  =  M  1z 


(E-2) 


Note  that  the  step  force  generates  both  homogeneous  and  a 
particular  parts  in  the  response.  The  homogeneous  part  arises 
from  the  initial  accelerations  and  the  particular  part  arises  from 
the  force  term  F(t)in  the  equation  of  motion.  For  simplicity,  the 
particular  part  of  the  step  response  will  not  be  introduced  until 
after  the  homogeneous  solution  is  determined.  The  particular 
solution  and  its  derivatives  of  order  up  to  and  including  2+a-/3 
are  initially  zero  and  have  no  effect  on  the  form  of  the 
homogeneous  part. 

Based  on  eqn  31  the  initial  values  of  the  homogeneous 
parts  of  the  basis  functions  should  satisfy 


r  £ 
6 


M*az 

0 


V  0 


W0) 

■  y3(0)  • 

y2(0) 
yx(0) 


(E-3) 


Here  M  z  is  the  prescribed  initial  accelerations  and  is  the 
s  t 

2n  +  1  column  vector  from  the  bottom  of  the  array.  To 
determine  the  homogeneous  part  one  needs  to  solve  for 
y j ( 0 ) »  j  “  1.2,3,  m-N.  The  orthonormal  transformation  matrix 
for  the  expanded  equations  of  motion,  eqn  12,  has  the  property 
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where  [M  J  is  the  pseudo  mass  matrix  associated  with  the  expanded 
equations  of  motion,  eqn  12.  Hence  premultiplying  both  sides  of 

7  ic  - 1 

equation  E-3  by  [$]  [M  ],  which  is  [4>]  ,  produces 
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Due  to  the  systematic  structure  of  [M  ]  and  the  eigenvectors 

associated  with  the  expanded  equations  of  motion  shown  in  Appendix 

“1  St 

B,  and  the  location  of  M  z  in  the  array  (the  2n+l —  column 
vector);  it  follows  that  E-5  reduces  to 
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from  which  one  concludes 
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where  q  is  defined  in  eqn  10.  Using  eqn  B-2,  this  expression  for 
the  initial  values  of  the  basis  functions  becomes 
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The  initial  values,  y  (0) ,  are  now  uniquely  determined; 
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however,  certain  conditions  resulting  from  eqn  E-8  will  be  useful 
in  simplifying  the  expressions  for  the  homogeneous  part.  Using 
eqn  E-8  to  substitute  for  the  initial  values  in  eqn  E-3  yields 
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and  as  a  result  one  sees  that 


m-N  -  , 


p  =  1,2,3,  ,m-l  p  h  2n  +  1  (E-10) 


and  again  using  eqn  B-2,  eqn  E-10  becomes 
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These  conditions  will  be  applied  to  the  homogeneous  part  of  the 
response . 

Based  on  eqn  31,  the  homogeneous  part  of  the  step  response  is 


m  •  N 

u(t)  =  I  .  y. (0) 
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(E- 12) 


Using  eqn  E-8  to  insert  the  initial  values  and  presenting  the 
Mittag-Leff ler  functions  in  terms  of  their  series  representations 
produces  the  homogeneous  part  of  the  step  response,  u  (t)  . 

— h 


u  (t) 

— h 

Z 


m  •  N 


J 


l  ±  b  ( -a^)q" 1  ^  z  l 
-1  J  3  3  i- 0 


(-aJ)V* 

r(i  +  *0) 


(E- 13) 


61 


Interchanging  the  order  of  summation  and  re-arranging  terms 
results  in 
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This  form  of  the  response  allows  one  to  apply,  in  a 
straightforward  manner,  the  conditions  given  in  eqn  E-ll,  which 
removes  several  of  the  lower  order  terms  in  the  time  series.  The 
resulting  form  of  the  homogeneous  part  of  the  step  response  is 
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Having  applied  all  the  initial  conditions  and  constructed  the 
homogeneous  part  of  step  response,  one  takes  the  first  derivative 
of  eqn  E-15.  Recall  that  (2n+q) P  is  equal  to  2+a. 
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Using  more  compact  notation  eqn  E-16  becomes 
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This  equation  is  the  expression  for  the  part  of  the  impulse 
response  dependent  on  the  initially  induced  velocities  (or 
momentum).  Recall  however,  that  particular  part  of  the  step 
response,  based  on  eqns  46  and  D-4, 
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is  not  yet  included.  Taking  its  first  derivative 
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and  combining  it  with  eqn  E-17  yields  eqn  47. 


APPENDIX  F  -  A  Class  of  Generalized  Functions 

Oldham  and  Spanier  (17:122)  present  the  one -half  order 
analogue  of  the  exponential  function.  This  function  is  predicated 
on  the  definition  of  fractional  differentiation  given  in  eqn  1, 
and  the  function  exhibits  singular  behavior  at  t  -  0 .  Note  that 
the  formulation  of  the  fractional  order  initial  value  problem  led 
to  the  exclusive  adoption  of  the  modified  definition  of  fractional 

A 

order  differentiation,  D[y(t)]  defined  in  eqn  23.  As  a  result 
one  is  motivated  to  construct  a  set  of  generalized  functions  based 
on  the  modified  definition  of  fractional  order  differentiation. 
As  expected  the  singular  behavior  does  not  appear  in  the 
derivatives  of  the  modified  analogues.  The  absences  of  the 
singular  behavior  lead  to  a  modified  set  of  generalized  functions 
having  properties  remarkably  similar  to  their  common  counterparts. 
Moreover,  this  set  of  modified  functions  is  precisely  those 
functions  found  in  the  solutions  to  the  modified  basis  equations 
which  are  the  foundation  of  the  fractional  order  initial  value 
problem . 

This  class  of  functions  is  built  around  the  special 
Mittag-Leffler  function 
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and  the  modified  definition  of  fractional  differentiation. 
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The  special  Mittag-Leffler  function  has  the  property 
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which  is  of  course  analogous  to 
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-r—  e  =  -  ae 
dt 


(F-4) 


In  fact  eqn  F-3  becomes  eqn  F-4  when  (3  is  set  to  one.  (13:327) 
Consequently,  for  our  purposes  the  special  Mittag-Leffler 
function,  eqn  F-l,  is  taken  to  be  the  generalized  fractional  order 
exponential  function.  The  property  given  in  F-3  may  be  proven  by 
taking  the  fi  order  derivative,  eqn  F-2,  of  each  term  in  the  series 
appearing  in  eqn  F-l. 

Having  established  che  fractional  order  exponential  function, 
the  definitions  of  the  fractional  order  sine  and  cosine  functions 
take  the  form 
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These  two  functions  have  derivatives  very  similar  to  the  regular 
sine  and  cosine  functions 
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Note  that  setting  /3  to  one  in  the  above  four  equations  yields  the 
exponential  representation  of  the  normal  sine  and  cosine  functions 
as  well  as  the  properties  of  their  derivatives.  Using  eqns  F-5 
and  F-6  one  can  straightforwardly  demonstrate  that 
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which  is  the  generalized  form  of  Euler's  formula.  This 

relationship  also  leads  to 
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which  for  /3  set  to  one  produces 

2  2 
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(F- 11) 


The  definitions  for  the  fractional  order  hyperbolic  sine  and 
cosine  functions  are 
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where  their  derivatives  are  seen  to  be 


sinh^  (at) 


cosh. 


(at) 


a^cosh 


=  a^sinh. 


,((at>'9) 


(F- 14) 


(at) 


t^ 


(F- 15) 


66 


These  functions  also  satisfy 


E^f±(at)^  =  cosh^  (at)^  ±  sinh^|(at)^j  (F-16) 

and. 

(at/  —(at/  =  cosh^(at)^]  -  sinh^(at)^  (F-17) 

All  of  these  definitions  and  relationships  for  the  generalized 
hyperbolic  functions  are  seen  to  reduce  to  traditional  definitions 
and  relationships  when  /3  is  one. 

The  remarkable  similarity  between  the  behavior  of  these 
generalized  functions  and  their  ordinary  counterparts,  not 
exhibited  by  Oldham's  generalized  functions,  is  directly 
attributable  to  the  absence  of  the  singular  behavior  in  Oldham's 
functions . 
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